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We study the finite-temperature transition to the m — 1/2 magnetization plateau in a model 
of interacting S = 1/2 spins with longer range interactions and strong exchange anisotropy on 
the geometrically frustrated Shastry-Sutherland lattice. This model was shown to capture the 
qualitative features of the field- induced magnetization plateaus in the rare-earth tetraboride, TmB4. 
Our results show that the transition to the plateau state occurs via two successive transitions with 
the two-dimensional Ising universality class, when the quantum exchange interactions are finite, 
whereas a single phase transition takes place in the purely Ising limit. To better understand these 
behaviors, we perform Monte Carlo simulations of the classical generalized four-state chiral clock 
model and compare the phase diagrams of the two models. Finally, we estimate a parameter set 
that can explain the magnetization curves observed in TmE>4. The magnetic properties and critical 
behavior of the finite-temperature transition to the m = 1/2 plateau state are also discussed. 
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I. INTRODUCTION 

Geometrically frustrated interactions and quantum 
fluctuation inhibit the stabilization of classical orderings. 
They sometimes become a trigger for the emergence of 
several exotic orders such as, quantum spin ice on the 
pyrochlore lattice[l|, spin liquid state on the honeycomb 
lattice , and multiple magnetization plateaus Q . There- 
fore, quantum spin systems with frustrated interactions 
have attracted great interest from both theoretical and 
experimental approaches. 

The 5 = 1/2 antiferromagnetic Heisenberg spin model 
on the Shastry-Sutherland lattice (SSL) Q is one of such 
systems. The Hamiltonian can be expressed by the near- 
est neighbor (intradimcr) J and the next nearest neigh- 
bor (interdimer) J' couplings: 

N.N. N.N.N. 

Experimentally, SrCu(B03)2[l, @ has been studied ex- 
tensively for its realization of the model. In this com- 
pound, CUBO3 layers stack along the c-axis direction 
and each magnetic layer consists of Cu 2+ ions carrying 
5 = 1/2 spins arranged in an orthogonal dimer structure 
that is topologically equivalent to the SSL. From sev- 
eral experimental observations, it was confirmed that the 
field dependence of the magnetization exhibits multiple 
magnetization plateaus. In these magnetization plateau 
states, Wigner crystals of spin-triplet dimersQ are real- 
ized reflecting the strong competition between the kinetic 
energy gain and mutual repulsion of the dimers. 

Theoretically, this model has been studied in great 
detail jl] and several states, such as the plaquette singlet 
state at zero field [p, [Io| and a spin supersolid state in the 
magnetic fields [ll[ were predicted. In a recent paper by 



Sebastian et al. [3j , the possibility of fractional magnetic 
plateaus was discussed in analogy to the quantum Hall 
effect. 

Fractional magnetic plateaus have recently been dis- 
covered at low temperatures in rare-earth tetraborides 
RB4 [R is a rare-earth element] [T34I1] • The magnetic 
moment carrying R 3+ ions in these compounds are ar- 
ranged in a SSL in the afr-planc. In TmB 4 , an extended 
magnetization plateau at m = 1/2 was confirmed for 
H cl ~ 1.9[T]< H < H c2 ~ 3.6JT] when a magnetic field 
is applied along the c-axis[l4lll7j. Here m is the normal- 
ized value by the saturation magnetization, m = m z /m s . 
In contrast to SrCu(B03)2 a strong anisotropy along the 
c-axis is expected owing to the crystal fields. From spe- 
cific heat measurements [13], it has been suggested that 
the degeneracy of the J=6 multiplet of Tm 3+ is lifted 
- the lowest energy state for a single ion is the non- 
Kramers doublet with J z =±6 and there exists a large 
energy gap to the first excited doublet. By restricting 
the local Hilbert space to the lowest energy doublet, the 
low-energy magnetic properties of the material can be de- 
scribed by a S=l/2 XXZ model with Ising-like exchange 
anisotropy and a ferromagnetic transverse coupling as 
discussed in the next section. 

The magnetization curves for the effective Hamilto- 
nian have already been calculated fl9l - [2~fl ] . The effective 
Hamiltonian is the S — 1/2 Ising-like XXZ model on the 
SSL and it is described by 

H = ]£(JS i -S 3 ) Ax + J2 ( j ' s >- s i)a ± 

N.N. N.N.N. 

- gf x B Hj2Si Z , (2) 

i 

where A_l denotes the Ising anisotropy and (JijSi • 
Sj) A± = -IJi^iSfSJ + SfS?) + JijSfS?. In the 
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Ising limit Aj_ = 0, it has been established by sev- 
eral approaches - such as Monte Carlo simulations [l9| 
and tensor renormalization-group analysis [2lj - that only 
m = 1/3 plateau is stabilized. The presence of the 
m = 1/2 plateau has been confirmed when quantum 
spin fluctuation is included fl9L [20I] . and the ground-state 
phase diagram for Aj_ > has been calculated. However, 
even for finite Aj_, the m = 1/3 plateau phase extends 
over a wider range of applied fields than the m = 1/2 
plateau phase. Since the to = 1/3 plateau is not ob- 
served in TmB 4 , the above Hamiltonian is insufficient to 
explain the experimental observations in T111B4. 



In a previous letter [22j, we argued that ferromagnetic 
J4 and antiferromagnetic J 3 couplings (see Fig. [T|) are 
necessary to explain the stabilization of an extended 
m = 1/2 plateau in the absence of the m = 1/3 plateau. 
We also investigated the finite-temperature phase transi- 
tion to the m = 1/2 plateau state. The results of finite- 
size scaling analysis indicated that a two-step second- 
order transition takes place - the difference between two 
critical temperatures is within 0.5% of J. The universal- 
ity class at both critical points is explained by the critical 
exponents of the two-dimensional Ising model. 

For the finite-temperature transition, a nontrivial 
question remains unanswered. In the m = 1/2 plateau 
phase, the lowest energy state is four-fold degenerate. 
Therefore, it is naively expected that the universality 
class should be the same as that of the four-state Potts 
model. As wc discuss in this paper, the critical behavior 
can indeed be explained by the four-state Potts univer- 
sality in the Ising limit, while the quantum spin model 
shows a two-step transition with both transitions belong- 
ing to the two-dimensional Ising universality class. Hence 
the phase diagrams for the thermal phase transitions for 
the two are different although both models possess the 
same symmetry. The low energy behavior of the m = 1/2 
plateau is shown to be described by the generalized four- 
state chiral clock model. As far as we know, the finite- 
temperature transition of the generalized chiral four-state 
clock model has not been studied precisely. Therefore, 
clarification of the finite-temperature transition to the 
m = 1/2 plateau phase is also valuable from the view 
point of statistical mechanics. In the present paper, we 
discuss the properties of the model introduced in our pre- 
vious letter in greater details, and clarify the nature of 
the phase transitions to the m = 1/2 plateau state. 

This paper is organized as follows. In section II, we re- 
produce the derivation of an effective Hamiltonian that 
captures the unique magnetization features in TmB4. An 
approach from the Ising limit of the effective Hamilto- 
nian shows that J4 couplings are important to stabilize 
the to = 1/2 plateau state. In section III, we focus on the 
finite-temperature properties of the effective model and 
re-investigate the universality class of the transition to 
the to = 1/2 plateau phase by extensive quantum Monte 
Carlo simulations. Wc show that the critical properties 
of the quantum model is different from those in the Ising 
limit. In section IV, a classical model, the generalized 




FIG. 1: Effective model on the SSL with diagonal coupling 
J, nearest-neighbor coupling J', and additional couplings J3 
and J4. 



chiral four-state clock model, is proposed in order to un- 
derstand the critical behavior of the original model. From 
Monte Carlo simulations, we find that the topology of 
the phase diagram for the classical model is the same as 
that of the original model. In comparison with the phase 
diagram for the classical model, we discuss the effects 
of quantum exchange on the finite-temperature transi- 
tion. In section V, we discuss the magnetic properties 
of TmB4 and predict the critical properties of the finite- 
temperature transition to the to = 1/2 plateau phase. 
Sec. VI is devoted to a summary of the results. 



II. EFFECTIVE MODEL 

In this section, we derive the effective model for the 
low energy part of magnetic properties of TmB4 again. 
As discussed in the introduction, the ground state of the 
single Tm 3+ ions is a J z = ±6 doublet with a large energy 
gap between the first excited doublet. At low tempera- 
tures (compared to this energy gap), one can ignore the 
contributions from the higher multiplets and the low en- 
ergy properties of the material are well described by an 
effective two-state model consisting of the J z = ±6 dou- 
blet. The fact that the interactions between the moments 
are derived from itinerant electrons enables us to expect 
isotropic type (Heisenberg type) interactions. There- 
fore, we start from the Hamiltonian, H' e ^ = H a + H' , 
where H = -D^J?) 2 , H' = V ,, ( ;,,.J,.I ,. and 
< \Gij\ <C D, where G%j is the coupling constant be- 
tween two moments on the site i and j, and D is the 
easy-axis anisotropy from the crystal fields. By treating 
H a as the biggest term, we find that there are four de- 
generate states, namely \J*,J*) = |6,6), | — 6,6), |6, — 6), 
and I —6,-6), in the lowest energy level. We apply per- 
turbation theory for H' e ^ and calculate the matrix el- 
ements among these four states. The Hamiltonian for 
arbitrary moment pairs can be described by the S =1/2 
Ising-likc XXZ Hamiltonian, hij = (JijSi ■ Sj) A , where 

(JijSi ■ S 3 ) A± = -\Jy\Aj. (SfSJ + SfS]) + JySfS*. 
Most importantly, the matrix elements for the transverse 
coupling ((6, —6\H'\ — 6, 6) and (—6, 6|i/'|6, —6)) are pro- 
portional to —{Gij/(—D)} 2J . Consequently, if the mag- 
netic moment J is even, the transverse coupling always 
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becomes ferromagnetic. (Note that the longitudinal com- 
ponent of the interaction becomes ferromagnetic or anti- 
ferromagnetic depending on the sign of Cry.) 

The interaction Gy is expected to be the RKKY type, 
because this compound is a metal. Therefore, the ef- 
fect of the long-range interactions is significant when the 
magnetic properties of TmB4 are considered. The prin- 
cipal interactions necessary to reproduce the dominant 
features of magnetization in TmB4 - in addition to the 
conventional SSL model with J and J' couplings - are 
the J 3 and J4 couplings, shown in Fig. [1] (f). There is 
another coupling with shorter range than that of J3 and 
J4 - the next-nearest neighbor coupling J 2 that is orthog- 
onal to J in the plaquettes with the diagonal interactions 
in the original SSL model. However, as we show in the 
following analysis, J 2 is not efficient in stabilizing the 
m = 1/2 plateau state, because it stabilizes the m = 1/3 
plateau at the same time. 

The Hamiltonian considered here is described by 

h = ( JS * ■ s , W + E ( J ' s * • s ; W 

<»>i>" <*.i>"' 
- gUBH^Si', (3) 

i 

where (ij), (ij)' , (ij)" , and (ij)'" denote sums over all 
pairs on the bonds with the J, J', J3, and J4 couplings, 
respectively. The positive (negative) sign of each cou- 
pling denotes antiferromagnetic (ferromagnetic) interac- 
tion. The original SSL model interactions are always 
assumed to be antiferromagnetic, i.e., J > and J' > 0. 
In the following, we set J as the unit of energy and ex- 
press all the parameters of the model in units of J. We 
studied the above model on square lattices of the form 
L x L with periodic boundary conditions. 

When A_l -C 1, the magnetic properties of the Hamil- 
tonian (J3|) are qualitatively explained by considering the 
Ising limit (Aj_ = 0). At J3 = J4 = 0, there are 
two candidates of the ground states at m — m z /m s = 
2N- 1 J2 i S z i = 0, where N = L d and d — 2. For 
J' I J < 0.5, each spin pair located on the diagonal bond 
J forms an antiferromagnetic classical dimer (ACD) 
(Fig. [5] (a)), but macroscopic degeneracy remains due 
to cancellations among interactions on the coupling J'. 
The local energy of such state can be estimated as 
Eacd = —J/8. In the other limit J'/ J 1, the sys- 
tem approaches the antiferromagnetic Ising model on the 
square lattice, and the Neel state becomes the ground 
state (Fig. [2] (b)). The local energy of the Neel state 
is calculated as E^eei = J/8 — J'/2. As is expected, 
the energy levels of the two states cross at J' / J = 0.5 
and T = 0, and the system shows a phase transition at 
the point. The local energies of the 1/2 and 1/3 plateau 
state are estimated in the same manner. The spin con- 
figurations in Fig. [21(c) and (d) are realized in the 1/3 
and 1/2 plateau phase, respectively. (These configura- 



tions were also confirmed from the snapshot of spin con- 
figuration in the Monte Carlo simulations in our previ- 
ous study [22I].) From the spin configurations, we obtain 
E 1/3 = - J/24 - J'/6 - H/6 and E 1/2 = -H/A, where 
£7i/3 and Ei/ 2 denote the local energy of the 1/3 and 
1/2 plateau state. Since the energy of the fully polarized 
state is given by Ep = J/8 + J'/2 — H/2, the magnetiza- 
tion curve shows a jump from the m = 1/3 to m = 1 at 
H c = J/2 + 2 J'. We show the field dependence of local 
energy for these five states in Fig. [3] The energies of the 
m = 1/2 and m = 1/3 plateaus, and the fully polarized 
state are always degenerate at the saturation field H c . 
Consequently, the degeneracy between these three states 
should be lifted when the additional couplings J 3 and J4 
are included. By estimating the energy gains due to J3 
and J4, we find that the ferromagnetic J4 coupling is one 
of the most efficient ways to stabilize the 1/2 plateau [27| . 
This analysis further demonstrates that the coupling J2, 
which is perpendicular to the diagonal coupling J, does 
not stabilize the m = 1/2 plateau, because the number of 
antiparallcl spin pairs on the J 2 couplings equals that of 
the parallel ones. Hence this coupling was not included 
in the effective Hamiltonian. 




FIG. 2: Conventional spin configurations in (a) an antiferro- 
magnetic classical dimer (ACD) state, (b) the Neel state, (c) 
1/3 plateau state, and (d) 1/2 plateau state. Solid and open 
circles indicate down and up spins respectively, (e) Schis- 
matic spin configuration of C2-symmetric phase (see text). 
Solid (open) ellipses mean spin pairs having the total magne- 
tization S z ~ (S z ~ 1). 



III. QUANTUM MONTE CARLO 
SIMULATIONS 



The qualitative discussion of the previous section is 
supplemented by large scale numerical simulation of the 
underlying model. We developed a variant of the stan- 
dard Stochastic-Series-Expansion quantum-Monte-Carlo 
method based on the modified directed-loop al gori thm 
to treat the longer-range interactions efficiently. [23l . [2f| 
Using this new algorithm, simulations of the Hamiltonian 
((21) were performed on square lattices L x L. Note that 
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FIG. 3: (color online) Field dependence of the local energy 
for the Ising limit. 



the negative sign problem does not appear in the quan- 
tum Monte Carlo computations for the Hamiltonian ([3]), 
because the transverse coupling is ferromagnetic as dis- 
cussed in the previous section. 

Fig. 2] shows the magnetization curves up to L = 40, 
for different values of the ratio J'/ J, and how the curves 
evolve upon the introduction of ferro- and antiferro- 
magnetic J4. In the computations, we fixed the J4 cou- 
pling constant at \J 4 /J\ = (l/\/2) 3 x J'/ J ~ 0.106 and 
the Ising anisotropy Aj_ at Aj_ = 0.2. 

At J3 = J4 = 0, two magnetization plateaus appear at 
to = 1/2 and 1/3, consistent with previous studies [1- 
[2~j| . When the longitudinal component of J4 coupling is 
ferromagnetic, the system shows a strong tendency to- 
wards the stabilizing an extended plateau at m = 1/2. 
Snapshots of the Monte Carlo simulations in the to = 1/2 
plateau phase confirm the realization of spin configura- 
tion shown in Fig. [2] (d). As long as there is the finite 
ferromagnetic J4 coupling, the to = 1/3 plateau region 
shrinks, for all values of the ratio J'/ J. This follows 
naturally from the analysis of the Ising spin model dis- 
cussed above. In the Ising limit, the field extent of the 
to = 1/2 plateau expands by an amount proportional to 
4| J4I/J, while that of the m = 1/3 plateau contracts by 
an amount proportional to 6| J^/J. 

When the longitudinal component of J4 coupling is 
antifcrromagnetic both to = 1/2 and 1/3 plateaus disap- 
pear. For antiferromagnctic J4, no conclusive evidence 
of any plateaus has been obtained in our calculations 
except for J'/ J = 0.5. A feature appears in the magneti- 
zation curve around to = 1/4 for J'/ J = 0.5. Fig. [5] (a) 
shows the expansion of the magnetization curve around 
H ~ 0.75 for J' I J = 0.5. For 0.60 < H/J < 0.74, 
the magnetization curve becomes flat at m = 2/9. Sig- 
nificantly, the to = 2/9 plateau was also obtained for 
the parent Shastry- Sutherland model and observed in 
SrCu(B03)2.[I[ To identify the spin structure, we cal- 
culated the spin-spin and the bond-spin correlation func- 
tions. (As for the definition of bond spins, please see eq. 
(01) The results are shown in Fig. [S] (b) and (c), respec- 
tively. In this plateau region, the 3^x3^ structure ac- 
companying 90 rotational symmetry breaking is clearly 
stabilized. This is the same symmetry breaking as the 



to = 1/2 plateau state. This plateau state has a charac- 
teristic feature: the periodicity of the to = 2/9 plateau 
state is longer than the distance of J4 couplings. Further- 
more, the periodicity of this plateau is the same as that 
discussed in ref. [2J], but the precise spin configuration 
is a little bit different from their results shown in figure 
5 in ref. (2~H because the magnitude of the moments on 
the qlaquette without the diagonal bonds is lager than 
that on the diagonal bonds having total S z = 1. 
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FIG. 4: Magnetization curves when J 4 /J = (1/V2) 3 J'/J 
(color online) and Ax = 0.2. Circles indicate the results at 
J4 = case. Triangles and squares are the results of the fer- 
romagnetic and antiferromagnetic longitudinal J 4 couplings, 
respectively. All symbols accompany error bars, which are 
smaller than the symbol size (here and the following figures). 

In the following, we focus on the finite-temperature 
transition to the to = 1/2 plateau state. The snapshot 
of the spin configuration of the to = 1/2 plateau phase 
leads us to expect the 90 (C4) rotation symmetry break- 
ing around the center of the plaqucttc without diagonal 
coupling J. In context of the bare spin language, the 
universality class of the finite-temperature transition is 
expected to be the four-state Potts universality class be- 
cause the lowest energy state in the to = 1/2 plateau is 
four-fold degenerate. However, the other scenario is also 
possible for A^ 7^ because the group C4 has a subgroup 
C2 ■ By introducing the 90 lattice rotation "C4" around 
the center, we can express the symmetry group in the 
paramagnetic phase as C4 = {e, C4, c\, C4 1 }. If the sys- 
tem exhibits a two-step phase transition, the symmetry 
breaks down to {e, c|} at the higher critical temperature 
T c i, and the remaining symmetry breaks down to the 
trivial group {e} at the lower critical temperature T C 2. 

To investigate the symmetry breaking at the critical 
point, we introduce two order parameters. The order 
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FIG. 5: (color online) (a) Magnetization curves when J'/J = 
0.5 and the antiferromagnetic J4/J = (1/V2) 3 J'/J. (b) 
and (c) are results for spin-spin and bond-spin correlations 
at fcflT/J = 0.02. Colored (open) circles correspond to the 
positive (negative) amplitude. The area of circles is propor- 
tional to the amplitude. Error bars are smaller than the sym- 
bol sizes, (d) Schematic spin configuration in the m = 2/9 
plateau state. Solid circles denote up spins. Yellow circles 
surrounded by up spins are occupied by up or down spins 
with almost equal probability. 

parameter of 90 rotation symmetry breaking can be ex- 
pressed by 

B* = {^\Y,(- 1 ) nrd)B ( r *)\), ( 4 ) 

where Yd is position of the diagonal bond, Bljd) is a 
product of longitudinal component of two spins on the 
bond Yd, and /(r^) takes ±1 depending on the position of 
the diagonal coupling (Fig. [6] (a)). That of 180 rotation 
symmetry breaking is also given by 

< = (^£{Sf(r)-S?(r)}>, 

r 

m % = <^£{S£(r)-Sf(r)}>, (5) 

r 

where the suffixes of longitudinal spin operators repre- 
sent the site indexes shown in Fig. [6](b), r is a positional 
vector for plaquettc. For T > T cl , the system is in the 
paramagnetic phase, and B st = and m x = are sat- 
isfied. For T cl > T > T c2 , the system retains the 180° 



(C2) rotation symmetry, and then B st 7^ and m x = 0. 
Below T C 2, the lowest energy state is characterized by 
B st + and m c x ^ 0. 

(a) (b) 




FIG. 6: (a) A bond spin operator. The bond spin operator 
B is defined by inner product of pair spins on each diagonal 
bond (bold line) and B 3t is the staggered magnetization of B. 
The bold black and gray lines indicate the sign of the factor 
f(rd)- (b) Definition of sublattice. 

First, we show the results of the finite-temperature 
transition in the Ising limit (Aj_ = 0). The Monte Carlo 
simulations up to L = 120 were performed at J'/J = 0.3 

and J 4 /J = -(l/V^) 3 x J'/J 0.106. Fig. (a) 

displays the temperature dependence of the correlation 
ratio of the bare spins at (x, y) = (L x /2, 0) and (L x /3, 0) 
for different system sizes. The curves cross at the crit- 
ical temperature T c ~ 0.09050(10), reflecting the size 
invariance of the correlation ratio at the critical tem- 
perature. The spin fluctuations freeze at this temper- 
ature. Fig. [7] (b) and (c) show the temperature de- 
pendence of the Binder ratios Rb for B st and Rs for 
mj, respectively. The system size dependence of both 
Binder ratios disappears at the same critical tempera- 
ture. The crossing of the curves for Rb and Rs also in- 
dicates that the C4 symmetry breaks down to the trivial 
group at T c ~ 0.09050(10). The obtained results suggest 
that the transition belongs to the universality class of the 
four-state Potts model. This is confirmed by finite-size 
scaling analysis. For the fourth order cumulant of B st 
and m x , we assume the scaling forms Ur = Fi^tL 1 /") 
and Us = FitL 1 /"), where F is a scaling function and 
Ur and Us denote the fourth order cumulant of B st 
and m c x , respectively. For the static structure factor of 
the bond spins, SbL 213 ^ — FftL 1 /") is assumed, where 
Sb = Ylr B(rd) exp [iTr(r x + r y )}. In the analysis, the 
critical exponents and the critical temperature are fixed 
at v = 2/3, = 1/12126] and T c = 0.09050(10), The 
results are shown in Fig. [8l The excellent data collapse 
confirms that the transition in the Ising limit belongs to 
the universality class of the four state Potts model. 

Next, we show the results of the quantum spin model 
with A_l = 0.25. We performed computations up to L = 
84. Fig. [5] shows the temperature dependence of the 
same order parameters as those defined in the above. The 
obtained results clearly indicate that the transition is a 
two-step process. 

The curves for R B cross at k B T cl /J = 0.0815(3). This 
means that Neel order of the bond spins B is realized 
for T < T c i. Since a quotient group C4/C2, which is 
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FIG. 7: (color online) Temperature dependence of the cor- 
relation ratio at J4/J ~ —0.106 in the Ising limit case, (a) 
Correlation ratio of the bare spin correlation, (b) The Binder 
ratio Rb of the bond-spin operator R at . (c) The Binder ratio 
Rs of m c x . 



isomorphic to the group C 2 , breaks down at T cl , it is 
naively expected that the transition belongs to the two- 
dimensional Ising model universality class. Accordingly 
we performed finite-size scaling analysis assuming the 
scaling form Ur = F(tL v l v ') and v = 1. We empha- 
size here that the data collapse shown in Fig. [5] (c) and 
(d) was obtained without any adjustable parameter, con- 
firming that the transition at T c ± indeed belongs to the 
two-dimensional Ising universality class. 

Fig. [9] (b) shows the temperature dependence of the 
Binder ratio Rs for m%, and it provides an evidence of 
the other phase transition. The curves for Rs for differ- 
ent system sizes intersect at fcsT c2 /J = 0.0580(5) and 
the difference T c \ — T c2 is approximately 0.03J, slightly 
larger than the previously obtained value (22|. The uni- 
versality class of the phase transition at T c2 should be 
the same as that of the two-dimensional Ising model be- 
cause the remaining symmetry C 2 is also isomorphic to 
Z 2 . This is confirmed by the finite-size scaling results 
presented in Fig. |9] (e) where we have used the same 
scaling function and critical exponent as for Rb, viz. 
U s = F(tL^ v ) and v = 1. Consequently, we conclude 
that there exists an intermediate phase with (^-rotation 
symmetry phase that can be characterized by the bond 
Neel order accompanying the internal antiferromagnetic 




FIG. 8: (color online) Finite-size scaling analysis for Ur, Sb, 
and Us in the Ising limit case, (a) and (c) are the results for 
the fourth-order cumulant Ur and Us- (b) Scaling analysis 
for the static structure factor Sb of bond spins. 



bare-spin fluctuation. The schematic spin configuration 
in this ^-rotation symmetry phase is shown in Fig. m 
(e). 

If the antiferromagnetic fluctuation on the bond J is 
relevant to the stabilization of the intermediate phase, 
the temperature range of the C2-rotation symmetric 
phase expands upon increasing Aj_. In Fig. [TU1 Aj_ de- 
pendence of the critical temperatures is presented for 
fixed J'/J = 0.3, H/J = 1 and J 4 /J = -(l/%/2) 3 x 
J 1 1 J ~ —0.106. The results show that the critical tem- 
perature associated with the breaking of the 180 rota- 
tional symmetry shifts to the lower temperature as Aj_ 
increases. This suggests that the dimerization of antipar- 
allel spin on the diagonal bonds plays a key role in the 
transition. As discussed in the followings, the phase dia- 
gram is qualitatively understood from a comparison with 
that of the generalized chiral four-state clock model. 



IV. COMPARISON TO A CLASSICAL MODEL 

In the previous section, we have investigated the finite- 
temperature phase transition to the m = 1/2 plateau 
phase for the classical and quantum spin models. Since 
the symmetry of both models is same and the finite- 
temperature phase transition is focused on, it is expected 
that the schematic phase diagram should be the same. 
However, the results obtained are different from this ex- 
pectation. In order to understand the discrepancy, we 
consider the generalized chiral four-state clock model. 
As we show below, this simplified classical model cap- 
tures the characteristics of phase diagram for the original 
model. 

We begin by introducing eight-spin unit clusters (wind- 
mill clusters) in the m = 1/2 plateau state as shown in 
Fig. Elb). In this phase, a down spin occupies one of 
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FIG. 9: (color online) Finite-temperature transition to the 
1/2 plateau phase, (a) and (b) are temperature dependence 
of the Binder ratio of Rb and Rs, respectively, (c) and (d) are 
the results of the finite-size scaling analysis for fourth-order 
cumulant Ur of staggered magnetization B 3 t and the static 
structure factor Sb of bond spins, (e) and (f) correspond to 
(c) and (d) for m%. In (f), the static structure factor of m% is 
calculated from Ss = ^ r m £(r) exp[— i{irr x + 7rr y )] 



four sites on a plaquette without the diagonal bond and 
up spins occupy the remaining three sites. If a down spin 
occupies the site labeled by '2' in Fig. [6jb) , the nearest 
down spin prefers to locate at the third-neighbor distance 
from it as long as the ferromagnetic J4 coupling. This 
indicates that the down spins in the m = 1/2 plateau 
state occupies the same edge of diagonal bond connect- 
ing each plaquette. Since the ferromagnetic J4 coupling is 
essential for stabilizing the m = 1/2 plateau state (as dis- 
cussed in Sec. II) , it is reasonable to describe the m = 1/2 
plateau state by an arrangement of such clusters. This 
description becomes exact in the limit \ J±/ J\ 3> 1 and 
Aj_ = 0. 

In the Ising limit Aj_ = 0, the spin configuration in 
the m = 1/2 plateau state can be expressed by the ar- 
rangement of the windmill clusters shown in Fig. 1111 (a). 
There are four kinds of clusters corresponding to the 90 °- 
rotation symmetry breaking in the m — 1/2 plateau. We 
assign clock spins pointing at an angle to each clusters. 
Here we ignore clusters having m 7^ 1/2 for simplicity. 
Thus we obtain the simplified classical Hamiltonian that 
can be regarded as the generalized four-state clock model 
on a square lattice with the nearest and next-nearest 
neighbor interactions, 

Hsim = Hn.n. + Hn.n.n 
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0.3 

Phase diagram at J'/ J = 0.3 and J4/J = 
~ —0.106. I.M. in the phase diagram means the 
C2-rotational symmetric phase. Schismatic spin configura- 
tions in the spin ordered phase and I.M. are shown in Fig. [2] 
(d) and (e), respectively. Open and solid circles indicate the 
second order transition with the two-dimensional Ising uni- 
versality, and an open square is that with the four-state Potts 
universality. An open triangle is the Kosterslitz-Thouless 
transition. All lines are guide to the eyes. 



%N.N.N = 7^ J'a 



JiJyVi + lJ + lj 



E 



where 6ij denotes an angle of clock spins and takes 
0,7r/2,7r, or 37r/2. The value of 9ij denotes the posi- 
tion of a down spin in the windmill clusters (see Fig. 
ITT1 (a) ) . J x and J y are interactions between two clock 
spins, and their values are listed in Table I. The interac- 
tions, J x and J y , are calculated from the sums of cou- 
pling energy between spins located in different clusters. 
When = 0, the values, A, B, C, and D, in J x and 
J y (see Table I) are evaluated as A=4J4, B=2J', C=0, 
D=4J', respectively. The resulting classical Hamiltonian 
retains C4 rotational symmetry conditionally: the sys- 
tem requires a lattice rotation when a global angle shift 
of all clock spins is executed. Such requirement arises 
from the geometric characteristics of the SSL. Since the 
symmetry of the Hamiltonian "Hn.n. is lower than that 
of "Hn.n.n., we focus on the critical properties of Hn.n. 
in the followings. 

The simplified Hamiltonian "Hn.n. is identical to 
the conventional four-state Potts model in the limit 
B=C=D=0. Therefore, it is trivial that the finite- 
temperature phase transition at this point is a sin- 
gle second-order transition belonging to the four-state 
Potts universality. However, the universality class be- 
comes nontrivial, when B^D^O. Performing Monte 
Carlo simulations, we numerically investigated the finite- 
temperature transition of the classical Hamiltonian 
"Hn.n. with D=2B and A=-4 up to L = 72. We sum- 
marize the obtained critical temperatures in Fig. Q~2] 
(a). These critical temperatures were estimated from the 
finite-size scaling analysis for the Binder ratio and corre- 
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TABLE I: Coupling constants of the simplified classical 
Hamiltonian. Note that (A, B, C, D) = (4J 4 ,2J',0, 4J') corre- 
sponds to the parameter set we have treated in the section III. 
Upper (lower) side tables are matrices of the nearest neighbor 
(next-nearest neighbor) interaction. 
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FIG. 11: (a) Windmill clusters and clock spins in the simpli- 
fied classical model (color online). A cluster is constructed by 
eight spins. Red dashed lines drawing windmills correspond 
to the clusters and arrows in the windmills denote angles of 
clock spins that represent states of the clusters, (b) General- 
ized four-state clock model on a square lattice. Solid squares 
correspond to each clusters described in (a). Solid (dashed) 
lines denote the nearest (next-nearest) neighbor interaction 
J x and J y (X and X). 

In the following, we focus on the results when C/A=0 
and B/A are varied (along the y-axis in Fig. [121(b) and 
(c)). For C/A=0, the system undergoes a single second- 
order transition when the coupling ratio satisfies B/A > 
—0.075. To identify the universality class, we performed 
finite-size scaling analysis and present the results in Fig. 
1131 (The scaling forms assumed here are the same as 
those shown in the previous section.) In the analysis, 
we estimated T c from the crossing of the curves for the 
correlation ratio, r coan $(L/2)/T cosn g(L/4), for different 
system sizes, where r cos „g(r) = (cosnOo.o ■ co$,n6 r ^ r ) and 
n = 1,2. For B/A > —0.07, excellent data collapse is 
obtained with v = 2/3 and j3 = 1/12. This strongly 
suggests that the transition belongs to universality class 
of the four-state Potts model[26j]. For B/A < -0.075, 
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FIG. 12: (color online) (a) Critical temperatures of the gener- 
alized four-state clock model when we fix at B/A=-0.05 and 
C/A=0. The higher and lower critical temperatures were cal- 
culated by the finite-size scaling analysis for the Binder ratio 
of (cos 6) and (cos 29) , respectively. Red and blue symbols 
correspond to the cross sections on red solid and dashed lines 
drawn in (b) and (c), respectively, (b) and (c) are expected 
phase diagrams of the present clock model. In (b) and (c), an 
area painted by green indicates the region where the two-step 
phase transition with two-dimensional Ising universality takes 
place. A blue quarter circle in (c) means the area where the 
system undergoes a single phase transition. 



the critical exponent increases from v = 2/3 to v = 1 
with the fixed ratio j3/v ~ 1/8 as B/A decreases. For 
B/A < —0.5, we confirm a clear evidence of two-step 
phase transition with v = 1 at both critical temperatures. 
Therefore, the critical properties for B/A < —0.5 belong 
to the two-dimensional Ising universality class. 

Next we discuss the results when the parameters are 
varied along the dashed lines in Fig. Q2] (b) and (c), 
where B/A=-0.05 and < C/A < 1. For C/A > 0.25, 
the system clearly undergoes a two-step phase transition 
with both transitions belonging to the two-dimensional 
Ising universality class. Fig. [represents the results of the 
finite-size scaling analysis at B/A=-0.05 and C/A=0.5. 
When the two-step phase transition takes place, the 90 ° 
rotation symmetry breaks at the higher critical temper- 
ature and the 180 ° rotation symmetry survives in the 
intermediate phase. The lower critical temperature de- 
creases as the difference between A and C decreases. For 
C/A > 1, the system seems to undergo a single phase 
transition and the 180 ° rotation symmetry survives for 
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any finite temperature. 
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FIG. 13: (color online) (al) - (a4) ((bl) 
results of finite size scaling analysis for the fourth-order 
cumulant of (cos 29) (the static structure factor of cos (9) 
when C/A=0. In these analysis, we estimated all critical 
temperatures from crossing points of the correlation ratios 
<r cosfl (L/2))/<r cose (L/4)> and (r cos 2fl (L/2))/(r cos 2fl (L/4)>. 
The data collapse in (al) and (bl) ((a4) and (b4)) was ob- 
tained by using the fixed critical exponents v — 2/3 and 
/3 = 1/12 (y = 1 and /3 = 1/8). The critical exponents in 
(a2), (b2), (a3), and (b3) were estimated to minimize the dis- 
crepancies among data. 



The obtained results allow us to consider two scenar- 
ios of the phase diagram for the present four-state clock 
model. One is that the system shows a single phase tran- 
sition only at B=C=0 and the two-step phase transi- 
tion takes place for B/A < and C/A > as shown 
in Fig. [T^] (b). The other is that there exists a re- 
gion having the four-state Potts universality class around 
B/A = C/A = 0, denoted by the area in Fig. QH (c). 
The critical exponents may show a crossover behavior at 
the boundary of the area. Further studies of the criti- 
cal behavior of the extended four-state clock model are 
currently underway. 

We compare the phase diagrams shown in Fig. [TUland 
Fig. [T2l here. The topology of the phase diagram for the 
simplified classical model is the same as that of the Ising- 
like XXZ model on the SSL. In both models, the region 
of the intermediate phase, where only the 90 ° rotational 
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FIG. 14: (color online) Finite size scaling analysis for fourth 

order cumulant of the order parameters (cos 29) (a) and 
(cos 6») (b) at B/A=-0.1 and C/A=0.5. (c) and (d) are the 
same analysis for the structure factor of each order parameter. 
In these analysis, each critical temperature is estimated from 
across of the correlation ratios (2 cos #l/2,o)/( cos ^L/i,o) and 
{cos9 L /2,o)/(cos8 L /4 i0 ). All results are obtained by using the 
fixed critical exponents v = 1 and f3 = 1/8. 



symmetry is broken, shrinks rapidly as the system ap- 
proaches the four-state-Potts universality point. When 
we consider the quantum effects in the original Hamil- 
tonian ([5]), fluctuation of antifcrromagnctic spin pairs 
on the diagonal bond is enhanced when increases. 
Then, the sublattice magnetization of antiparallel spin 
pair closes to zero. In the case, where the antiparallel 
spin pair becomes a classical or quantum-triplet dimcr 
with S z = 0, the clock spins pointing towards 9 = and 
7r (or 7r/2 and 3n/2) direction can not be distinguished 
from each other. Such dimcrization effect is captured as 
a reduction in the difference between A and C. Conse- 
quently, the lower transition temperature T c i associated 
with the 180 rotational symmetry breaking shifts to- 
ward T = as Aj_ increases. 



By considering the parameter set associated with the 
Ising limit of the original Hamiltonian, we can estimate 
a parameter set of 'Hn.n. as Asa 4J4, Bss 2J', C~ 0, 

and Dw 4 J'. For J'/ J = 0.3 and J 4 /J 0.106, Un.n. 

shows a single phase transition with the critical exponent 
v ~ 1, in clear disagreement with that of the original 
model. This is because we ignored the effect of the next- 
nearest neighbor couplings and the other states excluded 
from the four states treated here. A more quantitative 
estimation of the parameter set (A,B,C,D) is required 
by adding such effects for a quantitative mapping of the 
original Hamiltonian over the entire parameter range, but 
this is beyond the scope of this paper. 
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V. TmB 4 

Finally, we discuss the magnetization properties of 
T1T1B4. The to = 1/2 plateau state is observed for 
H cX ~1.9[T]< H <H c2 ~3.6[T] at 2[K][H|. The coupling 
ratio J'/J ~ 1 and a strong Ising anisotropy Aj_ << 1 
are expected from the cr ysta l structure of TmB4 and spe- 
cific heat measurements [18j. These experimental results 
help us estimate the other parameters for the present 
model. We obtain J 3 /J ~ 0.1182 and J 4 /J < -0.25 at 
J'/J = 1 from the local energy estimation in the Ising 
limit. The magnetization curves at fixed J3/J = 0.1182, 
Ji/J = -0.251, J'/J = 1 and k B T/J = 0.15 are shown 
in Fig. rT5](a). It is found that the magnetization jumps 
appear at H c \/J ~ 1.4 and H C 2/J ~ 2.65 not only in 
the Ising limit but also for A > 0. From these criti- 
cal fields, we roughly estimate {H C 2 — H c i)/H c i ~ 0.9. 
This value is in good agreement with the experimental 
value. Therefore our estimation for the coupling ratio 
is quantitatively consistent with the critical fields of the 
magnetization jumps observed in T111B4. 

The to ~ 1/8 plateau observed in the experiments 
[l|| Qj} appears to depend on the history of the system; 
to ~ 1/8 was observed around H ~ 1-8[T], when the 
fields were decreased from the saturation fields, whereas 
the magnetization remained vanishingly small around 
H ~ 1.8[T] during the upsweep. In other words, a hys- 
teresis loop was observed around H ~ 1.8[T]. The inset of 
Fig. [15] (a) is the magnetization curves at ksT/J = 0.05 
when we perform field sweeps. We started with the Neel 
state (the fully-polarized state) and increased (decreased) 
the applied field from H/J = (H/J = 3). A hysteresis 
behavior similar to the experiments can be observed in 
short-time simulations. Although the Monte Carlo dy- 
namics is not directly comparable to the real dynamics, 
this similarity is suggestive. The magnetization curves in 
the inset suggest a presence of shoulder around to = 1/8 
for 0.9 < H/J < 1.2. In this field region, we have 
confirmed the existence of a mixed state comprised of 
the Neel order and domain walls constructed of fully- 
polarized spin chains from the snapshot of spin config- 
uration. The period of the domain walls is fluctuating 
and the average of the period changes continuously as 
the field decreases. Thus, the shoulder around m = 1/8 
seems to be a transient process in successive plateaus 
having the magnetization < to < 1/2. Since the mag- 
netization value of the shoulder shifts to the lower mag- 
netization value as the system size increases, it seems to 
be smeared out in the thermodynamic limit. More quan- 
titative discussion of the to =1/8 shoulder is desirable via 
comparison with the other experimental observations. 

We comment the finite-temperature transition to the 
m = 1/2 plateau phase. Using the above parameters, we 
have studied it for the Ising spin model and the quan- 
tum spin model numerically. Fig. [15] (b) and (c) are the 
results of the scaling analysis, when the quantum spin 
model is considered. The results indicate that a single 
second-order transition takes place and the critical ex- 



0.75 



H/J=2.125 0>) 
k B T/J=0.233(3) . 
v=2/3 




0.25 



FIG. 15: (color online) (a) Magnetization curves at J'/J — 1, 
J3/J = 0.1182 Ja/J = -0.251. Red solid (open) symbols 
represent the results of the Ising limit case (the quantum spin 
model case) when they have been obtained by the field sweeps 
from H/J = 3 and the slow cooling. Blue ones are the results 
of the field sweeps from H/J — 0. Inset of (a) is the results 
of short time Monte Carlo calculations when the field sweep 
is performed, (b) and (c) show a finite-size scaling analysis 
for the fourth-order cumulant of Rb and Rs, respectively. 



ponent equals v ~ 2/3. The obtained results indicate 
that the critical property of the finite-temperature tran- 
sition to the to = 1/2 plateau phase can be explained by 
that of the four-state Potts model. This critical behavior 
was also confirmed in the Ising limit case. Therefore, if 
our estimation for the coupling constants is correct, we 
expect the experimentally observed critical exponents in 
TmB4 to belong to the four-state Potts universality class. 

The to =1/2 plateau has been observed in the other 
SSL compound ErLu[l3l|. The magnetic moments of 
ErB4 are derived from Er 3+ ions and the magnitude 
equals J = 15/2. In this compound, the strong Ising 
anisotropy perpendicular to the SSL plains has been sug- 
gested. When we focus on the magnetic properties at a 
very low temperature, the effective Hamiltonian of the 
magnetic properties can be also described by the present 
model in the Ising limit. Therefore we believe that our re- 
sults help understand the magnetic order of the to = 1/2 
plateau observed in ErB4. 



VI. SUMMARY 

In summary, we have studied the magnetic proper- 
ties of a model of interacting spins with Ising-like ex- 
change anisotropy and longer range interactions on the 
SSL that captures the low-temperature magnetic prop- 
erties of the rare-earth tetraboride, T111B4. We have fo- 
cused on the finite-temperature transition to the to = 1/2 
plateau phase and investigated the universality class of 
the transition. In the Ising limit, the system shows a sin- 
gle second-order transition. The critical behavior is well 
explained by that of the four-state Potts model. When 
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there exists quantum exchange interactions, it has been 
confirmed that there is the two-step phase transition with 
an intervening intermediate phase. Both the transitions 
are belonging to the two-dimensional Ising universality 
class. From the finite-size scaling analysis, we have as- 
certained that the intermediate phase can be character- 
ized by the " 180 "-rotation symmetric" state retaining 
7r-rotation symmetry on individual triplets with m z = 0. 
Since the symmetry of the quantum spin model is the 
same as that of the classical spin version, it is naively ex- 
pected that the critical behavior of both models will also 
be the same. However, the obtained results have indi- 
cated that the phase diagrams for the finite-temperature 
transition are different. To understand the critical behav- 
ior of both the models, we have proposed the simplified 
classical model, namely the generalized four-state chiral 
clock model. By performing the Monte Carlo computa- 
tions for the simplified classical model, we have found 
that the universality class at the critical temperatures 
and the topology of the phase diagram are in agreement 
with those of the original models on the SSL. Finally, 
we have studied the magnetic properties of TmB4. From 
the Ising limit analysis, we have estimated the parame- 
ters that can explain (qualitatively) the magnetization 



curves observed in experiments. From the short-time 
Monte Carlo simulations, we have suggested that the 
m ~ 1/8 plateau seems to be a metastable state. We 
have also investigated the finite-temperature transition 
to the m = 1/2 plateau state. If our estimation for the 
coupling constants is correct, measurements of the crit- 
ical exponents belonging to the universality class of the 
four-state Potts model are expected in experiment. 
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